'''
WARNING:
  Preserve Python3 compatibility
  print()
'''


'''
Jorg Stelling model from
A synthetic tunable mammalian oscillator
'''

import sys
import math
from numpy import *
from scipy import *
from bricks import *


BA_BISTABLE     = 0
BA_MONOSTABLE57 = 1
BA_MONOSTABLE97 = 2
BA_AMBIGOUS     = 3
BA_NOT_VALID    = -1


class BistableSwitch:
  '''
  ARGUMENTS
    TIME : simulation time (minutes)
  '''
  gDEBUG = True
  gTIME = 10000
  gG1 = 200.0
  gG2 = 200.0
  gG3 = 200.0
  gG4 = 200.0
  gG5 = 200.0
  gG6 = 200.0
  gpc_start = [];
  gpc_half_life = 0;
  gpc_duration = 0;
  grg_start = [];
  grg_half_life = 0;
  grg_duration = 0;
  gthresh = 1e2
  greq = False
  ghk = [];

  '''
  MODEL PARAMETERS


  G1: [57][59][pCMV][TAL97:KRAB][NEPTUN]
  G2: [59][97][pCMV][TAL57:KRAB][mCITRIN]
  G3: [PIP][pCMV][TAL97:KRAB]
  G4: [RIP][pmin][TAL57:KRAB]
  G5: [pCMV][PIP-KRAB]
  G6: [pCMV][Rapalog IP]
  '''

  gZERO = 1e-4
  


  def __init__(self, dbg, t, g1, g2, g3,g4,g5,g6, req, pcs, pchl, pcd, rgs, rghl, rgd, hk, thr):
    self.gDEBUG = dbg
    self.gTIME = t
    self.gG1 = g1
    self.gG2 = g2
    self.gG3 = g3
    self.gG4 = g4
    self.gG5 = g5
    self.gG6 = g6
    self.gpc_start = pcs;
    self.gpc_half_life = pchl;
    self.grg_start = rgs;
    self.grg_half_life = rghl;
    self.gthresh = thr;
    self.gpc_duration = pcd;
    self.grg_duration = rgd;
    self.greq = req;
    self.ghk = hk;



  #################
  # SIMULATION PART
  #################



  def plotResults(self,TIME, T97m, T97p, T57m, T57p, PBSm, PBSp, PBSa, RBSm, RBSp, RBSa, PC, RG, NEPm, NEPp, MCTm, MCTp, BIF97p, BIF97t, BIF57p, BIF57t):
    if sys.version_info[0] > 2:
      return

    from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
    from matplotlib.figure import Figure
    from matplotlib.dates import DateFormatter
    from random import random
    
    NONCESIZE = 1000;
    NONCE = int(random()*NONCESIZE)

    if self.gDEBUG:
      f0 = figure()
      f1 = figure()
      f2 = figure()
      f3 = figure()
      f4 = figure()
    else:
      f0 = Figure()
      f1 = Figure()
      f2 = Figure()
      f3 = Figure()
      f4 = Figure()
    
    s0=f0.add_subplot(111)
    s0.plot(TIME, RBSm, linestyle='dashed', label='RIP mRNA', color='blue')
    s0.plot(TIME, RBSp, label='RIP', color='blue')
    s0.plot(TIME, PBSm, linestyle='dashed', label='PIP mRNA', color='red')
    s0.plot(TIME, PBSp, label='PIP', color='red')
    s0.set_xlabel('Time [min]')
    s0.set_ylabel('Concentration [nM]')
    s0.legend(loc="upper right")
    f0.set_facecolor('white')
    lg = s0.legend(loc="upper right")
    lg.get_frame().set_alpha(0.5)    
        
    s1=f1.add_subplot(111)
    s1.plot(TIME, PBSa, label='PIP-active', color='blue')
    s1.plot(TIME, RBSa, label='RIP-active', color='red')
    s1.plot(TIME, PC, label='Pristinamycin', color='blue', linestyle='dotted')  
    s1.plot(TIME, RG, label='Rapalog', color='red', linestyle='dotted')
    s1.set_xlabel('Time [min]')
    s1.set_ylabel('Concentration [nM]')
    s1.legend(loc="upper right")
    f1.set_facecolor('white')   
    lg = s1.legend(loc="upper right")
    lg.get_frame().set_alpha(0.5)
     
    s2=f2.add_subplot(111)
    s2.plot(TIME, T97m, linestyle='dashed', label='TAL97:KRAB mRNA', color='blue')
    s2.plot(TIME, T97p, label='TAL97:KRAB', color='blue')
    s2.plot(TIME, T57m, linestyle='dashed', label='TAL57:KRAB mRNA', color='red')
    s2.plot(TIME, T57p, label='TAL57:KRAB', color='red')
    s2.set_xlabel('Time [min]')
    s2.set_ylabel('Concentration [nM]')
    s2.legend(loc="upper right")
    f2.set_facecolor('white')   
    lg = s2.legend(loc="upper right")
    lg.get_frame().set_alpha(0.5)
    
    s3=f3.add_subplot(111)
    s3.plot(TIME, T97m, linestyle='dashed', label='NEPTUN mRNA', color='blue')
    s3.plot(TIME, T97p, label='NEPTUN', color='blue')
    s3.plot(TIME, T57m, linestyle='dashed', label='mCITRIN mRNA', color='red')
    s3.plot(TIME, T57p, label='mCITRIN', color='red')
    s3.set_xlabel('Time [min]')
    s3.set_ylabel('Concentration [nM]')
    s3.legend(loc="upper right")
    f3.set_facecolor('white')   
    lg = s3.legend(loc="upper right")
    lg.get_frame().set_alpha(0.5)
    
    
    s4=f4.add_subplot(111)
    s4.plot(T57p, T97p, 'c.', label="Phase")
    s4.plot(BIF57p,BIF97p, 'rx', label="Stable state")
    s4.plot([self.gthresh for i in range(10000)],linspace(1e-2,1e4,10000), label="Threshold State #1", linestyle='dashed', color='gray');
    s4.plot(linspace(1e-2,1e4,10000), [self.gthresh for i in range(10000)], label="Threshold State #2", linestyle='dashed', color='gray');
    s4.loglog()
    s4.set_xlabel('mCITRIN')
    s4.set_ylabel('NEPTUN')
    s4.legend(loc="upper right")
    f4.set_facecolor('white')   
    lg = s4.legend(loc="upper right")
    lg.get_frame().set_alpha(0.5)
    
    
    
    if self.gDEBUG:
      path = './out/'+self.__class__.__name__
      paths = [path+'_s0_'+str(NONCE)+'.png', path+'_s1_'+str(NONCE)+'.png', path+'_s2_'+str(NONCE)+'.png', path+'_s3_'+str(NONCE)+'.png', path+'_s4_'+str(NONCE)+'.png']
      f0.savefig(paths[0], format ='png')
      f1.savefig(paths[1], format ='png')
      f2.savefig(paths[2], format ='png')
      f3.savefig(paths[3], format ='png')
      f4.savefig(paths[4], format ='png')
    
    else:
      path = '/var/www/imodel/img/'+self.__class__.__name__
      paths = [path+'_s0_'+str(NONCE)+'.png', path+'_s1_'+str(NONCE)+'.png', path+'_s2_'+str(NONCE)+'.png', path+'_s3_'+str(NONCE)+'.png', path+'_s4_'+str(NONCE)+'.png']
      canvas=FigureCanvas(f0)
      f0.savefig(paths[0], format ='png')
      canvas=FigureCanvas(f1)
      f1.savefig(paths[1], format ='png')
      canvas=FigureCanvas(f2)
      f2.savefig(paths[2], format ='png')
      canvas=FigureCanvas(f3)
      f3.savefig(paths[3], format ='png')
      canvas=FigureCanvas(f4)
      f4.savefig(paths[4], format ='png')
        
    return paths
    

  def bifurcationAnalysis(self,BIF97p, BIF97t, BIF57p, BIF57t, T97p, T57p):
    '''
    BA_BISTABLE     = 0
    BA_MONOSTABLE57 = 1
    BA_MONOSTABLE97 = 2
    BA_AMBIGOUS     = 3
    BA_NOT_VALID    = -1
    '''
    
    state57 = False
    state97 = False
    ret     = BA_NOT_VALID

    #Discover ambigous states
    for i in range(len(BIF57p)):
      a = (BIF57p[i] >= self.gthresh and BIF97p[i] < self.gthresh)
      b = (BIF57p[i] < self.gthresh and BIF97p[i] >= self.gthresh)
      if not (a or b):
        ret = BA_AMBIGOUS
        break;
      if BIF57p[i] >= self.gthresh:
        state57 = True
      if BIF97p[i] >= self.gthresh:
        state97 = True  

    if not ret == BA_AMBIGOUS:
      if state57 and state97:
        ret = BA_BISTABLE
      elif state57:
        ret = BA_MONOSTABLE57
      elif state97:
        ret = BA_MONOSTABLE97
      
    return ret
    
      
    
      
    
  def run(self):
    TIME = range(self.gTIME)
    #Gene dosage plasmid [plasmids/cell]
    #Relative number to cell 
    #TODO: to use with measurements transfer function
    #TODO: HEK 293
    #Batard 2001
    #*Assumption: uniform distribution
    #*Tunable
    G1 = self.gG1
    G2 = self.gG2
    G3 = self.gG3
    G4 = self.gG4
    G5 = self.gG5
    G6 = self.gG6
    



    #Experimental TAL Effector rates
    #RLA Measurements
    if self.gDEBUG:
      tal97 = Effector('TAL97:KRAB', 'dat/tal97krab.dat')
      tal57 = Effector('TAL57:KRAB', 'dat/tal57krab.dat')
    else:
      tal97 = Effector('TAL97:KRAB', '/home/root/cellulator-mirror/dat/tal97krab.dat')
      tal57 = Effector('TAL57:KRAB', '/home/root/cellulator-mirror/dat/tal57krab.dat')
    #tal97.plotTransferFunction()
    #tal57.plotTransferFunction()
    #show()
    #quit()
    
    r97 = (G3+G1)/(G2)
    r57 = (G2+G4)/(G1)
    
    
    t57max = self.gthresh*10;
    t97max = self.gthresh*10;
    hk97 = self.ghk[2]
    hk57 = self.ghk[3]
    if self.greq:
      q97 = 0.1
      q57 = 0.1
    else:
      q97 = tal97.getTransferRate2(r97)
      q57 = tal57.getTransferRate2(r57)
    #print q97
    #print q57
    k57 = 1.0/(1/(q57)-1);
    k97 = 1.0/(1/(q97)-1);


    #Inducer affinity
    #(testing)
    #TODO: look for source
    Apc = 1.0
    Arg = 1.0

    #Maximal transcription rates (mRNA/min)
    #apply pCMV, pMin
    #Apply protein length
    k1 = 30.0
    k2 = 30.0
    k3 = 30.0
    k4 = 30.0
    k5 = 30.0
    k6 = 30.0


    #Basal activity of inducible promoter
    #TODO: source, RLA, mammalian
    #TODO: Jan, Miha
    alfapc = 0.085 #8.5% represija
    alfarg = 0.025 #40x aktivacija

    #Binding affinity 
    #Literature data. PIP repressor.
    #Tigges.
    KPC = 3.0
    KRG = 3.0

    #Cooperativity
    #TODO: source
    npc = self.ghk[1]
    nrg = self.ghk[2]

    #Translation rates (protein/min)
    #TODO: source
    #TODO: protein length adjust
    t57 = 0.02
    t97 = 0.02
    tpc = 0.02
    trg = 0.02
    tnep = 0.02
    tmct = 0.02


    #Degradation rate (mRNA/min)
    #TODO: source
    kdm = 0.0173

    #Degradation rate (protein/min)
    #TODO: source
    kdp = 0.0058

    
    #Degradation rate (inducer/min)
    kdpc = 2/(exp(self.gpc_half_life))
    kdrg = 2/(exp(self.grg_half_life))
    
    #Concentration scaling factor for copies/cell
    #TODO: source
    fv = 1.85e-3*0.5

    #Temporary switch
    #TODO: sources
    r1 = k1*fv
    r2 = k2*fv
    r3 = k3*fv
    r4 = k4*fv
    r5 = k5*fv
    r6 = k6*fv

    '''
    OBSERVED SPECIES
    '''
    RBSm = array([0.0 for i in TIME]) #rapalog binding protein mRNA
    RBSp = array([0.0 for i in TIME]) #rapalog binding protein
    RBSa = array([0.0 for i in TIME]) #rapalog binding protein (active)

    PBSm = array([0.0 for i in TIME]) #pristinamicin binding protein mRNA
    PBSp = array([0.0 for i in TIME]) #pristinamicin binding protein
    PBSa = array([0.0 for i in TIME]) #pristinamicin binding protein (active)
    PBSad = array([0.0 for i in TIME]) #pristinamicin binding protein (active/draw)

    PC = array([0.0 for i in TIME]) #pristinamicin
    RG = array([0.0 for i in TIME]) #rapalog

    if self.gpc_half_life > 0:
      PC[self.gpc_start] = 1500
    else:
      l = len(self.gpc_start)
      for i in range(l):
        PC[self.gpc_start[i]:(self.gpc_start[i]+self.gpc_duration[i])] = 1500
        
    if self.grg_half_life > 0:
      RG[self.grg_start] = 1500
    else:
      l = len(self.grg_start)
      for i in range(l):
        RG[self.grg_start[i]:self.grg_start[i]+self.grg_duration[i]] = 1500
    
    T97m = array([0.0 for i in TIME])
    T97p = array([0.0 for i in TIME])
    T57m = array([0.0 for i in TIME])
    T57p = array([0.0 for i in TIME])
    
    NEPm = array([0.0 for i in TIME])
    NEPp = array([0.0 for i in TIME])
    MCTm = array([0.0 for i in TIME])
    MCTp = array([0.0 for i in TIME])
    
    BIF97p = []
    BIF97t = []
    BIF57p = []
    BIF57t = []

    '''
    SIMULATION
    '''

    
    for t in TIME:
      #Skip first and last steps
      if t == 0 or t==self.gTIME-1:
        continue

      
      #TODO: check with article
      dPBSm = G5 * r5 - kdm * PBSm[t-1]
      dPBSp = tpc * PBSm[t-1] - kdp * PBSp[t-1]
      dPBSa = PBSp[t-1] * (PC[t-1]/(Apc + PC[t-1]))
      dPBSad = PBSp[t-1] * (1 - PC[t-1]/(Apc + PC[t-1]))
      
      #TODO: check with article
      dRBSm = G6 * r6 - kdm * RBSm[t-1]
      dRBSp = trg * RBSm[t-1] - kdp * RBSp[t-1]
      dRBSa = RBSp[t-1] * (RG[t-1]/(Arg + RG[t-1]))

      #r3 = fv * k3 * (alfapc + (1-alfapc)*(1-PBSa[t-1]**npc/(KPC + PBSa[t-1]**npc)))
      alfapc = alfarg
      r3 = fv * k3 * (alfapc + (1-alfapc)*(PBSa[t-1]**npc/(KPC + PBSa[t-1]**npc)))
      r4 = fv * k4 * (alfarg + (1-alfarg)*(RBSa[t-1]**nrg/(KRG + RBSa[t-1]**nrg)))
    
      #TAL 97
      dT97m = G1 * r1 * 1/(1+(T57p[t-1])**hk57/(k57*t57max)**hk57);
      T97m[t] = T97m[t-1] + dT97m
      T97m[t] = int(T97m[t] > 0) * T97m[t]
      dT97m = G3*r3 - kdm * T97m[t-1]
      T97m[t] = T97m[t] + dT97m
      T97m[t] = int(T97m[t] > 0) * T97m[t]
      dT97p = t97 * T97m[t-1] - kdp * T97p[t-1]

      #TAL 57
      dT57m = G2 * r2 * 1/(1+(T97p[t-1])**hk97/(k97*t97max)**hk97);
      T57m[t] = T57m[t-1] + dT57m
      T57m[t] = int(T57m[t] > 0) * T57m[t]
      dT57m = G4*r4 - kdm * T57m[t-1]
      T57m[t] = T57m[t] + dT57m
      T57m[t] = int(T57m[t] > 0) * T57m[t]
      dT57p = t57 * T57m[t-1] - kdp * T57p[t-1]

      #NEPTUN
      dNEPm = G1 * r1 * 1/(1+(T57p[t-1])**hk57/(k57*t57max)**hk57)  - kdm * NEPm[t-1];
      dNEPp = tnep * NEPm[t-1] - kdp * NEPp[t-1]
      
      #MCITRIN
      dMCTm = G2 * r2 * 1/(1+(T97p[t-1])**hk97/(k97*t97max)**hk97) - kdm * MCTm[t-1];
      dMCTp = tmct * MCTm[t-1] - kdp * MCTp[t-1]

      #PRISTINAMICYN
      dPC = -PC[t-1] * kdpc
      
      #RAPALOG
      dRG = -RG[t-1] * kdrg
      
      
      #BIFURCATION
      #All possible phases with dT97p/dt = 0, dT57p/dt = 0
      if math.fabs(dT97p) < self.gZERO and math.fabs(dT57p) < self.gZERO :
        #Avoid initial state
        if T57p[t-1] > 0.0 or T97p[t-1] > 0.0:
          BIF97p.append(T97p[t-1])
          BIF57p.append(T57p[t-1])
        
      PBSm[t] = PBSm[t-1] + dPBSm
      PBSp[t] = PBSp[t-1] + dPBSp
      RBSm[t] = RBSm[t-1] + dRBSm
      RBSp[t] = RBSp[t-1] + dRBSp   
      PBSa[t] = PBSa[t] + dPBSa
      PBSad[t] = PBSad[t] + dPBSad
      RBSa[t] = RBSa[t] + dRBSa
      T97p[t] = T97p[t-1] + dT97p
      T57p[t] = T57p[t-1] + dT57p
      NEPm[t] = NEPm[t-1] + dNEPm
      NEPp[t] = NEPp[t-1] + dNEPp
      MCTm[t] = MCTm[t-1] + dMCTm
      MCTp[t] = MCTp[t-1] + dMCTp  
      
      if self.gpc_half_life > 0 and not t in self.gpc_start:
        PC[t] = PC[t-1] + dPC
      if self.grg_half_life > 0 and not t in self.grg_start:
        RG[t] = RG[t-1] + dRG

      #Preserve above zero
      PBSm[t] = int(PBSm[t] > 0) * PBSm[t] 
      PBSp[t] = int(PBSp[t] > 0) * PBSp[t]
      RBSm[t] = int(RBSm[t] > 0) * RBSm[t]
      RBSp[t] = int(RBSp[t] > 0) * RBSp[t]
      PBSa[t] = int(PBSa[t] > 0) * PBSa[t]
      PBSad[t] = int(PBSad[t] > 0) * PBSad[t]
      RBSa[t] = int(RBSa[t] > 0) * RBSa[t]
      T97p[t] = int(T97p[t] > 0) * T97p[t]
      T57p[t] = int(T57p[t] > 0) * T57p[t]
      NEPm[t] = int(NEPm[t] > 0) * NEPm[t]
      NEPp[t] = int(NEPp[t] > 0) * NEPp[t]
      MCTm[t] = int(MCTm[t] > 0) * MCTm[t]
      MCTp[t] = int(MCTp[t] > 0) * MCTp[t]
      PC[t] = int(PC[t] > 0) * PC[t]
      RG[t] = int(RG[t] > 0) * RG[t]
      
      T97p[t] += int(T97p[t] == 0) * float(1e-2)
      T57p[t] += int(T57p[t] == 0) * float(1e-2)

      
    if self.gDEBUG:
      #print self.bifurcationAnalysis(BIF97p, BIF97t, BIF57p, BIF57t, T97p, T57p)
      return self.bifurcationAnalysis(BIF97p, BIF97t, BIF57p, BIF57t, T97p, T57p)
      pass
    return self.plotResults(TIME, T97m, T97p, T57m, T57p, PBSm, PBSp, PBSad, RBSm, RBSp, RBSa, PC, RG, NEPm, NEPp, MCTm, MCTp, BIF97p, BIF97t, BIF57p, BIF57t)


      

